Designing arrays of Josephson junctions for specific static 

responses 

J.G. Caputo 1 ' 2 and L. Loukitch 1 

1) Laboratoire de Mathematiques, INS A de Rouen 
B.P. 8, 76131 Mont-Saint-Aignan cedex, France 
2) Laboratoire de Physique theorique et modelisation, 
Universite de Cergy-Pontoise and C.N.R.S., France 

February 6, 2008 



Abstract 

We consider the inverse problem of designing an array of superconducting Josephson junc- 
tions that has a given maximum static current pattern as function of the applied magnetic 
field. Such devices are used for magnetometry and as Terahertz oscillators. The model is a 
2D semilinear elliptic operator with Neuman boundary conditions so the direct problem is 
difficult to solve because of the multiplicity of solutions. For an array of small junctions in a 
passive region, the model can be reduced to a ID linear partial differential equation with Dirac 
distribution sine nonlinearities. For small junctions and a symmetric device, the maximum 
current is the absolute value of a cosine Fourier series whose coefficients (resp. frequencies) 
are proportional to the areas (resp. the positions) of the junctions. The inverse problem 
is solved by inverse cosine Fourier transform after choosing the area of the central junction. 
We show several examples using combinations of simple three junction circuits. These new 
devices could then be tailored to meet specific applications. 

1 Introduction 

The coupling of two Type I superconductors across a thin oxide layer is described by the two 
Josephson equations pQ, 

F = $ ^, / = sJ c shu». (1) 

where <j) is the phase difference between the two superconductors in units of <J> = h/2e the reduced 
flux quantum, V and / are respectively, the voltage and current across the layer, s is the contact 
surface and J c is the critical current density. The Josephson equations and Maxwell's equations 
imply the modulation of DC current by an external magnetic field in the static regime (SQUIDs) 
and the conversion of AC current in microwave radiation [2 [3] . In all these systems there is a 
characteristic length which reduces to the Josephson length Aj, the ratio of the electromagnetic flux 
to the quantum flux $o for standard junctions. The behavior of a Josephson junction depends on 
its size compared to Aj. In small junctions the phase will not vary much except for large magnetic 
fields. Long junctions on the contrary enable large variations of the phase accommodating so-called 
"fluxons" or sine-Gordon kinks where the phase varies by 2ir [2]. 

For many applications and in order to protect the junction, Josephson junctions are embedded 
in a so-called microstrip line which is the capacitor made by the overlap of the two superconducting 
layers. This is the "window geometry" where the phase difference satisfies an inhomogeneous 2D 
damped driven sine Gordon equation [3] resulting from Maxwell's equations and the Josephson 
constitutive relations |T]). For resonator applications this design allows to couple the junctions in an 
array to increase the output power and adapt impedance for coupling the device to a transmission 
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line. In addition one can select some desirable dynamic features like resonances [8 j and optimize 
the frequency response over a given band for wave mixing applications [5]. 

Parallel arrays of Josephson junctions can be used in the static regime as very fine magnetic 
field detectors. The maximal current I max which can cross the device (see Fig. [1} for a given 
magnetic field H, without any voltage (V = the static regime) defines the I m ax(H) curve. The 
behavior of arrays of identical and equidistant small Josephson junctions has been extensively 
studied [21 [3] . The problem of finding I max remains difficult to solve because of the multiplicity 
of solutions due to the sine nonlinearity and the Neuman boundary conditions. 

For fundamental reasons and applications it is interesting to work with non-uniform arrays 
where the junction sizes and their spacings can vary. In [H[9] we developed a continuous/discrete 
or long wave model where the phase variation is neglected in the junctions and where the couplings 
between junction and surrounding microstrip are correctly taken into account. In particular we 
consider the waves between the junctions that are completely neglected by the classical Resistive 
Shunted Junction (RSJ) lumped models Our approach allows to choose the distance between 
junctions and their area. In the same device we can model junctions with different areas and 
different current response, in particular 7r-junctions. This simple model allows to analyze in depth 
the statics of the device and this is not possible from the 2D original equations [9]. This long 
wave approximation can be generalized to 2D to explain the behavior of squids [13j . In addition 
we obtain an excellent agreement with the complex experimental I max (H) curves [10] using the 
very simple magnetic approximation introduced in [S]. 

For experimentalists, it is very useful to extract parameters of the array from the I m ax(H) 
curve. For example it gives informations on the quality of the junctions. Recent studies by Itzler 
and Tinkham examine how defects in the coupling affect this maximum current [141 115] . This is 
important because high T c superconductors can be described as Josephson junctions where the 
critical current density is a rapidly varying function of the position, due to grain boundaries. 
Fehrenbacher et al[16] calculated J max (ff) for such disordered long Josephson junctions and for 
a periodic array of defects. The expressions obtained are complicated so the inverse problem of 
determining junction parameters from the 7 max (-ff) curve is very difficult to solve for arrays or 
general current densities. However, when the simple magnetic approximation of the I max (H) holds, 
it allows to extract information on the sizes and positions of the junctions in an array assuming 
Imax(H) is a periodic and even function. This is the purpose of this article. In particular we will 
show how one can obtain a cosine profile and multi-cosine profile from a combination of simple 
3 junction arrays. We will indicate what parameters can be obtained from a general I m ax{H) 
profile. After presenting the general model in section 2, we introduce the magnetic approximation 
and give its properties in section 3. Section 4 discusses the inverse problem for a three junction 
array. In section 5 we design the device from a general I m ax(H) and conclude in section 6. 

2 The model 

The device we model (see Fig.QJ is a so-called microstrip cavity (grey area in Fig.rjJ) between two 
superconducting layers containing small regions (junctions) where the oxide layer is very thin (~ 
10 Angstrom) enabling Josephson coupling between the top and bottom superconductors. The 
dimensions of the microstrip are about 100 fim length and 20 /im width and the length and width 
of the junctions is about wj — 1 /im. In the static regime, the phase difference if between the 
top and bottom superconducting layers obeys the following semilinear elliptic partial differential 
equation [3] 

- Atp + g(x; y) sin <p = 0, (2) 

where g(x; y) = 1 in the Josephson junctions and outside and where we have neglected the differ- 
ence in surface inductance between the junction and passive region. This formulation guarantees 
the continuity of the normal gradient of (p, the electrical current on the junction interface. The 
space unit is the Josephson length Aj, the ratio of the flux formed with the critical current density 
and the surface inductance to the flux quantum $o- 
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Figure 1: The left panel shows the top view of a superconducting microstrip line containing 
three Josephson junctions, H, I and </> are respectively the applied magnetic field, current and 
the phase difference between the two superconducting layers. The phase difference <j> between the 
two superconducting layers satisfies — A<fr — in the linear part and — A</> + sin(</>) = in the 
Josephson junctions. The right panel shows the associated 2D domain of size I x w containing 
n = 3 junctions placed at the positions y = w/2 and x = Oi, i = 1, n. 



The boundary conditions representing an external current input / or an applied magnetic field 
H (along the y axis) are 
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where < v < 1 gives the type of current feed. The case v = 1 shown in Fig. [T] where the current 
is only applied to the long boundaries y = 0, w is called overlap feed while v — corresponds to 
the inline feed. 

We consider long and narrow strips containing a few small junctions of area wf placed on the 
line y = w/2 and centered on i = a„ i — l,n as shown in Fig. [1] Then we search tp in the form 
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where the first term takes care of the y boundary condition. For narrow strips w < it, only the 
first transverse mode needs to be taken into account [5] because the curvature of ip due to current 
remains small. Inserting ((4]) into ([2]) and projecting on the zero mode we obtain the following 
equation for tfio where the index has been dropped for simplicity 
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and where 7 = I/w and the boundary conditions 4>'(0) = H— {l—v) r y/2, and (f>'(l) = H+{l—v) r )/2. 
The factor Wi/w is exactly the "rescaling" of Xj(= 1) into A e // = > 1 due to the presence 

of the lateral passive region [7j. 

As the area of the junction is reduced, the total Josephson current is reduced and tends to 
zero. To describe small junctions where the phase variation can be neglected but which can carry 
a significant current, we introduce the following function gh 



Wj 



gii(x) = — , for a, — h < x < a,; + h, and gh(%) = elsewhere, 
where i = 1, ..n. In the limit ft-tOwe obtain our final delta function model [S] 
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where di = wf /w and the boundary conditions are 

0'(O) =H-(1- i/) 7 /2, 4>'{l) =H+(1- !/) 7 /2. 

Despite its crude character the delta function approximation is a good model for arrays with short 
junctions as long as the magnetic field is small compared to 1/wi where Wi is the size of the 
junctions ID] , It allows simple calculations and in depth analysis that are out of reach for the 2D 
full model. In addition when di < the model can describe so-called 7r-junctions. For these, the 
tunneling current is sin(</> + n) = — sin(0) instead of the usual sine term in the second Josephson 
equation |T]) . This new type of coupling occurs in some materials [111 112] and it is hoped to be 
incorporated in the design of arrays. It is then natural to associate negative di coefficients to n 
junctions in the device. 

We have the following properties. 

1. Integrating twice (0 shows that the solution <fi is continuous at the junctions x = a^, i = 
1, . . . n. 

2. Almost everywhere (in the mathematical sense), —<j)"(x) — vj/l, so that outside the junc- 
tions, <f> is a piece-wise second degree polynomial, 

4>{x) = -i^-x 2 + BiX + Ci , Va; e]a,;a i+ i[. (8) 



3. At each junction (x — flj), <f>' is not defined, but choosing ei > 0, and 62 > 0, we note that 
lim lim / 4>"{x)dx = \ (x)dx = [<f>'{x)} a >_ . 



Since the phase is continuous at the junction x = di, (pi = <p{ a i) we § e t 

[4>'{x)] a ji =d i sm{ ( j )i ). (9) 



4. Integrating |7]) over the whole domain, 

a n 
W\q = / 4>"dx = y^di sin(0i) - z^7 , 

i=l 

and taking into account the boundary conditions, we obtain 

a 

7 = y^^sin(^), (10) 

i=l 

which indicates the conservation of current. 

In [9j, we developed two ways to find the ^maxiH) curve for the device using these properties, see 
the Appendix " Piece- wise polynomial" for more details on the solution of the problem. The most 
useful property in [9] for this study is the magnetic approximation of the j m ax(H) curve. 



3 The magnetic approximation 

a+ 

Since [4>'\ i_ = di sin(0j) (remark[9|) and 7 < di, we notice that for small di, </> tends to the linear 

function 4>{x) = Hx + c. Starting from 4>{x) = Hx + c, it is simple to find the j m ax{H) curve. This 
is what we call the magnetic approximation. We generalize here what was done in [18] for arrays 
of equidistant junctions. We have shown in [9] that the j ma x{H) curve of ((7]) tends to it when 
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Figure 2: We plot "f m ax{H) and c max {H) for two four junction devices. In the left panel, we have 
a circuit with equidistant junctions of equal di = 1. For the right panel, the junctions are such 
that h = §, = |j ^3 = 1) with d% — 1.2, d 2 = 1.5, e?3 = 0.5, e?4 = 0.8. Notice that c max (H) 
varies in a complicated way for a non uniform device. 



di tends to zero. In experiments the di coefficients are small enough so that this approximation 
is justified and provides a quantitative estimate of the "f m ax(H) curve [lOj . For inhomogeneous 
arrays of many junctions, this curve is complex and even in this case the approximation is very 
good. 

Since 4>{x) = Hx + c then 7 = di sin(iJ<2i + c). To find the 'jmax(H) curve of the magnetic 
approximation, we take the derivative of 

/ n \ / n \ 

7 = ( di sin(iZai) J cos(c) + ( dj cos(Hai) J sin(c) = Acos(c) +Ssin(c), (11) 

with respect to c where we have isolated the factors A, B. The values of c such that dj/dc = 
are 

c max (H) = arctan f Sr 1 j^"^ ) , (12) 
and as we want a maximal (not only an extremal) current we obtain: 



2J sin (Ha t + c max (H)) 



(13) 



Now, we focus on the case where c max (H) is not defined. In this case, considering previous equation 
(fl^)) . we obtain Ym=i ^* sin(Hai) = = A. From Asin(c) = Bcos(c) we obtain cos(c) = or 
B = 0. Note that cos(c) = imply 7 = 0, in the other hand, B = 0, imply dj/dc — whatever 
the value of c. So, 7 is constant, and 7 ma: c = A = 0. We plot j m ax(H) and c max {H) in Fig.@, 
for a uniform device and for a non uniform one. In the second case, we have chosen the position 
of the junction to have a long period (H p = 12w). We notice that the length I of the device does 
not appear in eq. (|13p . 

In order to study the function (|13p . we start with a few definitions. 



U: We define the distance between consecutive junctions: 

H = a i+l ~ a i ■ 

Junction unit: We call a junction unit the set of distances between junctions. We denote it 

{h] h', ln-l} ■ 

We define the position of the junction unit as a±, the position of the first junction. 
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Figure 3: Curves j v 
d,2 = 1, d,2 — 2 and di 



versus H for three junction devices, from left to right di = 0, 
3. For all panels, a% = —1, a-i = 0, 03 = 1 and d\ — d-$ = 1. 



if,: For an n-junction device, Zf, = a„ — ai is junction unit length. 
Symmetric unit: We call a symmetric unit, an n-junction circuit such that 
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In the Appendix we prove the following three Propositions 

Invariance by translation (|7.1| : 7 ma3: (i?) does not depend on the position of the junction unit. 
Parity of the j max curve (|7.2[) : 7 maa: (i?) is an even function of H . 



Solution for a symmetric device (|7.3|) : For a symmetric junction unit such that a„ = —ax, 
c max (H) = ±tt/2 (this is not obvious from fig. (J2J))- 

In these propositions, we establish the most important result of this article. The ^ m ax curve 
for a symmetric junction unit can be calculated simply by centering the junction unit so that 
Cmax(H) — ±7r/2. More precisely, consider an n + 1 symmetric junction unit where n is even. We 
can always choose this by setting the d of the central junction {a\ + a n )/2 to 0. Then we shift the 
junction unit by a,( n +i)/2 so that the central junction is placed at x — 0. We relabel the junctions 
by setting i' = i — (n+l)/2. Then the central one is for i = 0, the first one to the right is i = 1, 
the first one to the left is i = — 1 ... so that the equation (|13|) becomes 
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(n+l)/2 

2 di cos(Hai) 



(14) 



where we omitted the primes. In the rest of the article we will consider the array to be symmetric. 



4 The direct problem for j ma x 
4.1 A device such that 'jrnax(H) = cos(if) 

In Fig. ([3]) we present from left to right the j m ax(H) for a SQUID (2 junctions), a uniform 3 
junction unit, a d\ = 1 = ^3, c?2 = 2 (termed 1-2-1) junction unit that is discussed in 3j and a 
1-3-1 junction unit. In all cases the junctions are equidistant. The first two panels, represent well 
known devices. 

Applying eq.fTTJJ to the following case di = 2, we obtain, 

lmax(H) = \d Q + di cos(iJai) = 2 + 2cos(i?ai) . (15) 
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Figure 4: Left panel: plot of j m ax{H) for a symmetric five equidistant junction device. The 
continuous line corresponds to different di while the dashed line is for equal di. See the text for 
the parameter values. The right panel shows the corresponding device. 

This is an exact cosine function shifted by a constant. With the last case d 2 = 3, we obtain 
lmax{H) = 3 + 2cos(-ffai). 

Comparing all the panels we understand the role of the central junction. We can have an exact 
representation of jmax for this type of circuit, if we imagine jmax as an absolute value of a simple 
cos function translated by the value do (which is equal to zero if there is no junction). Eq. (|f 4p 
shows that we can sum cosine functions, with a chosen amplitude and period. Remark that if 
d = —2 (it junction) then "f m ax{H) = 2 — 2cos(-£/ai). 

4.2 A multi-cosine 7 ma:r (i/) 

For arrays with more than two junctions, experimentalists can play on the set of distances U 
separating the junctions as well as on the strength di (proportional to the area) of each junction. 
We now show the influence of each set of parameters starting from the d^s. Fig. [4] presents 
on the left panel jmaxiH) for a symmetric set of 5 equidistant junctions a\ = 1, a 2 — 2. The 
dashed line corresponds to d» = 1, i = — 2 ... 2 giving a maximum current of 5. Here one sees the 
typical interference pattern between the main bumps. The small oscillations can be eliminated 
by choosing do = 1-82025, d\ = d_j = 1.25 and di = d- 2 = 0.3425 as seen from the continuous 
line on the left panel of Fig. [U This "pulse" profile could be very useful for specific applications 
because of the large region where 'ymax(H) = 0. The right panel of Fig. [4] presents what would 
be the device for this set of m and di. We chose a critical current density of 10 4 ^4 cm" 1 so that 
Xj ~ 10/im. We chose a transverse width w = lA^m. Assuming the area of the smallest junction 
to be 1/im, 2 , we get the scheme shown, where the central junction has an area 5.32 /im 2 . 

The second parameter that can be changed is the position a.- L of each junction in the array. 
As an example in Fig. [5] we show in the left panel the so-called "triangle" jmaxiH) obtained by 
setting ai = 1, a 2 — 3, d Q — 2.4888, d\ = 1.1234 and d 2 — 0.121. The dashed line presents 
lmax(H) for equal strengths. Changing the di's allows to eliminate the oscillations in the minima 
of "/max{H) and obtain an almost linear behavior. The right panel shows the arrangement of the 
junctions in the microstrip. We have chosen the same physical parameters as for Fig. 3J Now we 
can design all devices which have a j m ax(H) curve as a sum of di cos(aix), with di positive. We 
can notice that if a^/ai € K\Q, we can construct a non periodic 7 mal (]J). 

5 The inverse problem for a given 7 maa; (i^) = JgiH) 

We now show how to design an n + 1 junctions circuit (n is an even integer) to obtain a given 
7 g (H) curve. The formula (jT4j) can be used to solve this inverse problem using cosine Fourier 
transforms. To avoid ambiguities we assume a symmetric array and a positive and periodic j(H). 
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Figure 5: Left panel: plot of "f m ax(H) for a symmetric five junction device where the junctions 
are not equidistant, a\ = 1 and 02 = 3. The continuous line corresponds to different di while the 
dashed line is for equal di. The right panel shows the corresponding device. 



We have the following result. 

Proposition 5.1 (Solution of inverse problem for ^maxiH)) Assume aj g (H) even, periodic 
of period H p and strictly positive. The array is harmonic and the positions of the junctions are 
given by a 2 : = i2ir / H p where i is an integer. Their strengths di are given by 

d_j =d l = -^- { " j g (H) cos{Ha t )dH , Vi G {0, . . . , n/2) . (16) 
h p Jo 

This gives the positions and coefficients di of an array that will have a j m ax(H) that is the 
truncation to order n of the cosine Fourier series of j g (H ) . 

To gain insight into the problem let us first review the "pulse" example studied in the previous 
section. Assume j g (H) to be the 2n periodic extension of e~ aH where a is large enough. The 
coefficients di are given by 

d l = — e- aH2 cos(Hi) dH+— e -*(^-H) 2 cos ( jff j) dH = - e- aH ' 2 cos(Hi) dH. 
2tt Jo 2tt J tt J 

(17) 

These Fourier coefficients decay exponentially as expected [17] because 7 9 (-ff) is C°° over the 
interval [0, 2ir] and satisfies the boundary conditions. This means that i < 2 is enough to get a 
good approximation of 'jg(H). In fact Fig. 0] corresponds to j g (H) « 5e~ 6H and the formula 
(IT71) gives the values do = 1.82025, d\ = d-i = 1.25 and c?2 = d-2 = 0.3425 that were obtained 
in the previous section. The next coefficients d3 = 0.043, d^ = 0.0023 are very small and can be 
neglected. 

Let us now consider a square "f m ax{H) curve which could make a very fine magnetic detector 
because of its strong response over a given interval and zero response elsewhere. For that we 
assume the the square profile 

~f g (H) = 1 for nil - — ) < H < tt ( 1 - — ] and elsewhere, (18) 



and extend it periodically every 2tt — H p . To compute the parameters ai and di, we apply the 
previous result (see ea. ([T6| ') to obtain 

If 2 * , , fiirH\ 2 ( h 1 + h 2 \ fin h 2 - hA , . 

4 = -/ , g (H) C0S — ^ = -sin (i,^ cos y + ^ (19) 
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Figure 6: We compare the "f m ax{H) of the magnetic approximation to the exact solution of 
equation for v = (inline current feed. The parameters are given in the text. 



This gives the following values of di for h\ = hi = h 
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Note the decay in 1/i of di because J g (H) is only continuous. Another interesting fact is that 
some di are negative so that some junctions are it junctions. So we obtain an array of eleven 
junctions whose positions arc given above together with their strengths di (positive for a normal 
junction and negative for a it junction). 

In Fig. [S]we plot the magnetic approximation and the exact solution of equation Q for v = 
(this case is called inline current feed, see the section "Piece- wise polynomial" in the Appendix 
and for more details [9]). The values are / = 20, the junction unit is shifted by 10, the Josephson 
characteristic length is Xj = 5.6/im so that all di are multiplied by 0.035714285. 

We see that for this type of junction (about lfim 2 of area) inline current feed for ([7]) and 
magnetic approximation give close results. Differences appear when the maximal current is larger 
but the Gibbs phenomena is less important in the solution of the equation ([7]) than magnetic 
approximation. 



6 Conclusion 

Using a simple approximation, we introduced a method to design a symmetric array of Josephson 
junctions which has a specific "f m ax(H) static response. The sizes of the junctions are given by the 
coefficients di of the cosine Fourier transform of j m ax{H). Their position is a, = i2n/H p where 
Hp is the period of "f m ax(H). We use In + 1 junctions to obtain a curve formed with n Fourier 
coefficients. 

This work follows closely the article [Hj, where all the mathematical results were established, 
in particular the convergence of the solutions of the full problem ([7]) to the ones obtained in the 
magnetic approximation. There we show that the overlap current feed can cause a non even 
lmax{H) (see proposition "Magnetic shift" in [S]). 

If we are in the region of validity of our original model, ie the magnetic field is small and the 
distance between junctions is larger than their size, then we can design a device for a given static 
response. 
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7 Appendix 



7.1 Propositions 



Proposition 7.1 (Invariance by translation) The j ma x(H) curve obtained within the mag- 
netic approximation does not depend on the position of the junction unit. 

Proof. Let us assume two devices with the same junction unit {li; I2; l n } with the first 
junction placed respectively at x = a\ and x = a\ + c. We note 7m aa: (-ff) (respectively 7m OI (ff)) 
the jmaxiH) curve of the first (respectively the second) device. In the same way we note c x max {H) 
(respectively c^ aa ,(iJ)) the c max {H) function of the first (respectively the second) device. 



7 2 

/ ni ax 



(H) 



n 

E 

;.=i 



d.sin (H ai + Hc + c 2 max (H)) 



we note c 1 {H) = Hc + c\ lax (H'). As we do not know if c 1 (_ff) = o^^H), c x max {H) being the best 
value (if it exists) to obtain the maximal 7. Consequently, r ^ nax {H') < 7™ a:E (-ff)- 
On the other side, considering 



E 

i=l 



di sin (H{a t + c) + c^^H) - He) 



noting c 2 (H) = c^^H) — He and using the previous argument, we obtain 7m OX (-ff) > Imaxi^)- 
From the two previous inequalities, we obtain: J^n ax (H) = 7^ aa .(i?). □ 

Thus, the "f m ax(H) curve for the magnetic approximation depends only on the junction unit. 

Proposition 7.2 (Parity of the j max curve) For all devices, "f m ax(H) = ^maxi—H). 

Proof. Since sin and arctan are odd functions and cos is an even function then c ma x{—H) = 
(H) (see (pj). Finally, 



A-h) 



di sin (-Ha,i - c max (H)) 



rXH) 



Notice that c max is an odd function and jmax is an even function (see Fig. ([2]) ). 

Proposition 7.3 (Particular solution for symmetric device) For all symmetric units such 
that a n = ~a x , c max (H) = ±n/2. 

Proof. To see this, relabel the junctions so that the central one corresponds to i = 0, the 1st 
on the left to i = —1, the 1st on the right to i = +1.. Using the first proposition we can shift the 
junction unit so that ao — 0. Then the total current is 



1/2 



di sin(Ha,i + c) 

i=-n/2 

n/2 

do sin(c) + 2 di sm(Hai + c) 

i=l 
n/2 

sin(c) I do + 2 2_. di cos(Hai 

i=l 



Then 



= 7r/2 when do + 2^™fJ d, cos(Hai) > and c„ 



-7r/2 otherwise. Thus, 



(n+l)/2 

do + 2 di cos(Hai) 



(20) 
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7.2 Piece-wise polynomial 

Let <fi be a solution of ([7]) and <\>\ — <fi(ai). From remark ©, </> is a polynomial by parts. We 
define Pi+i(x) the second degree polynomial such that Pi + i(x) = 4>(x) for < x < aj+i. Using 
the left boundary condition we can specify 4> on [0; a\\: 

Pi(x) = (x 2 - a 2 ) + - i^ 7 ) (* - ax) + ^ (21) 

At the junctions ((9]) tells us that Vfc G {1; . . . ; n}, 

Pk+M - Pk(a k ) = d k $m{P k {a k )). (22) 

Considering that cj>" — —vj/l on each interval, the previous relation and the continuity of the 
phase at the junction, we can give a first expression for Pk+i, 

P k+ i(x) = ~j~ {x - a k ) 2 + [P k (a k ) + d k smP k (a k )} (x - a k ) + P k (a k ). (23) 

So, cj> is entirely determined by <px, 7 and H. 

The polynomials (|2ip and ([23)) establish existence and shape of <f> at junctions. Let see boundary 
conditions. The first, 

0'(Q) = i*(0) = H _(l_»/) 7 /2 
is true by construction; the second (for n junction circuit) is: 

P' n+x {l)=H+{l-u)l, (24) 

is true only for solutions of Eq. (0 . At H given, solutions of Eq. (|24|) define a relation between <p\ 
and 7. 

So, the maximal current solution depend on <pi and 7, and Eq. (|24[) is the constraint for search of 
lmax{H)- As solutions 4> are defined at 2tt almost (see equation (J7J), we can assume <j>i G [— 7r; 7t]. 
In other hand, Rem. ([TO]) teach us that J ma x € [0; J^. dj], so we take 7 € [0; J^. d,]. To solve this 
problem with Maple ©, we plot implicit function (the constraint) of two variables <f>\ and 7, with 
and ^ fixed. Let us note all variables: 

P'n+i | a=I 7, u,B)-H- ^7 = 0. (25) 

with ((^1,7) G [— 7r;7r] x [OjJ^dj]. Lastly the program search in exhaustive way, the biggest value 
of 7 of this implicit curve. Incrementing H 7 we obtain "f ma x(H). 

This method has the advantage to converge to the global maximal "/max [H]- 



12 



